#!/usr/bin/python3
# -*- coding: utf-8 -*-
# @Author  : lkm

"""
    GM(1,1)模型预测：北方某城市1986年-1992年道路交通噪声平均声级数据预测
"""

import numpy as np

x = np.array([71.1, 72.4, 72.4, 72.1, 71.4, 72.0, 71.6])
k = 2
### 计算级比
lamda = x[:len(x) - 1] / x[1:len(x)]
# 计算级比边界
left_boundary = np.exp(-2 / (len(x) + 1))
right_boundary = np.exp(2 / (len(x) + 1))

""" 判断级比 """
if lamda.min() >= left_boundary and lamda.max() <= right_boundary:
    print('通过级比检验，可用原始数据进行建模！')
    """ 对原始数据进行一次累加 """
    x_cumsum = x.cumsum()

    """ 构造数据矩阵B """
    B0 = (x_cumsum[:len(x_cumsum) - 1] + x_cumsum[1:]) / 2.0
    B0 = B0.reshape((len(B0), 1))
    B = np.append(-B0, np.ones_like(B0), axis=1)

    """ 构造数据向量Y """
    Y = x[1:].reshape((len(x) - 1, 1))

    """ 计算：发展系数a，灰色作用量b """
    [[a], [b]] = np.dot(np.dot(np.linalg.inv(np.dot(B.T, B)), B.T), Y)

    """ 建立模型：往后预测2步 """
    # if k == 0:
    #     result = ((x[0] - b/a) * np.exp(-a * k)) + (b/a)
    # else: # 输出n+1位置的预测值
    #     result = (((x[0] - b/a) * np.exp(-a * (len(x)+k)) + (b/a))) - (((x[0] - b/a) * np.exp(-a * (len(x)+k-1))) + (b/a))

    """ 求已知数据的预测值 """
    pred_y = [((x[0] - b / a) * np.exp(-a * 0)) + (b / a)]
    for i in range(1, len(x)):
        pred_y0 = (((x[0] - b / a) * np.exp(-a * i)) + (b / a)) - (((x[0] - b / a) * np.exp(-a * (i - 1))) + (b / a))
        pred_y.append(pred_y0)

    """ 计算残差、相对误差、级比偏差"""
    epsion = x - pred_y  # 残差
    delta = epsion / x  # 相对误差
    rho = 1 - (1 - 0.5 * a) / (1 + 0.5 * a) * lamda  # 级比偏差

    """ 预测第八个数据1993年的值 """
    y1993 = (((x[0] - b / a) * np.exp(-a * 7)) + (b / a)) - (((x[0] - b / a) * np.exp(-a * (7 - 1))) + (b / a))

else:
    print('没有通过级比检验！')